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Abstract. We present a technically simple treatment of self-trapping of Bose- 
Einstein condensates in double well traps based on intuitive semiclassical 
approximations. Our analysis finally leads to a convenient closed form 
approximation for the time-averaged population imbalance valid in both the mean- 
field case and in the case of finite particle numbers for short times. 
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1. Introduction 

Despite its simplicity an atomic Bose-Einstein condensate (BEC) in a double well trap 
is a quantum system showing a rich dynamical behaviour in different parameter ranges 
that can be quite accurately controlled in current experiments (see e. g. [J [21 [5| and 
references therein). Even in the mean- field limit of high particle numbers N where 
many-particle correlations can be neglected, different dynamical regimes are observed. 
The latter can be understood by mapping the system to a classical pendulum the 
dynamics of which can be analysed in the corresponding two-dimensional phase space 
(see e.g. [4]). For small and moderate interaction U between the particles a condensate 
that is initially localized in one of the two wells performs Rabi respectively Joscphson 
oscillations between the two wells. If the inter-particle interaction exceeds a critical 
value, the initially occupied well remains macroscopically occupied due to spontaneous 
symmetry breaking. This is referred to as running phase self-trapping effect. 

While for finite particle numbers the self-trapping effect is eventually destroyed by 
quantum correlations in the long-time limit, it can still be observed for a considerable 
time span depending on the number of particles in the system [SJ |B] . In this "short 
time regime" the main dynamical features of the system are well described by (semi- 
) classical phase space methods like the truncated Wigner approximation [7] or a 
similar ansatz based on the Husimi distribution |8l [9] where the respective quantum- 
mechanical phase space representations of an initial quantum state are sampled by an 
ensemble of trajectories that are propagated according to the corresponding Gross- 
Pitaevskii mean-field equation. Related aspects of the system studied in the literature 
include the semiclassical WKB-type quantization of its energy spectrum |10j and the 
quantum fluctuations of the cigenstates near the point of transition to self-trapping 

In this paper, following the example of |121 113] . the transition to self-trapping 
is described by means of a single scalar quantity, namely the time-average z{A) of 
the relative population imbalance between the two wells as a function of the scaled 
interaction parameter A = U {N —1) / J where J is the tunneling coefficient (cf. equation 
([T|) below). It will be shown that a useful closed form approximation for ^(A) in 
the aforementioned semiclassical parameter and time regimes can be derived in a 
technically simple manner by means of intuitive heuristic approximations. 

This article is organized as follows: In section [2] we discuss some selected 
properties of the system Hamiltonian, namely the Bose-Hubbard dimer, and its mean- 
field limit. In section[3]we use a heuristic ansatz to derive a closed form approximation 
for the time-averaged population imbalance as a function of the scaled interaction 
strength for the mean-field case. In section |3] many-particle effects are taken into 
account in a semiclassical manner by incorporating quantum fluctuations of the initial 
state into the mean-field description derived in the previous section which again leads 
to a closed form approximation for the time-averaged population imbalance. The main 
results are briefly summarized in section [5] [Appendix A| contains a brief review of the 
dynamics of the corresponding nonintcracting system. 
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2. The Hamiltonian 

In the two mode approximation the double well system is described by the Bose- 
Hubbard dimer Hamiltonian (of. e.g. [3]) 

H = (aloR + aj^aL^ + ^ ((4)^«L + (or)^"!!,) + CLa^aL + CRflRaR (1) 

where ol and or are the annihilation operators of a particle in the left and right well 
respectively. J > determines the rate of tunneling between adjacent lattice sites, 
U is the inter particle interaction parameter and el, er are the on-site energy terms. 
Here, we concentrate on repulsive interaction U > 0. 

For a sufficiently high number of particles N in the system, one can apply a 
Hartree-Fock mean-field approximation (see e.g. jl4j ) where the operators ul.r in 
the Hamiltonian can be replaced by the complex numbers ^/Ncl^r representing their 
respective coherent state expectation values. This leads to the "classical" Hamiltonian 

H/N = (c*CR + c* cl) + |(iV - 1) (|cl|^ + IcrI'') + + erI^rP • (2) 

The canonical equations of motion ihci = dH/dc'^, ihc\ = —dH/dci then lead to the 
coupled discrete Gross-Pitaevskii equations or nonlinear Schrodinger equations 

*/iCL=-^CR+(eL + C/(A^-l)|cLp)cL (3) 

ihcn = -^CL + (er + U{N - 1)|cr|2) cr . (4) 

which describe the mean-field dynamics of the on-site amplitudes. Particle number 
conservation yields 

Introducing the amplitude-phase representation cr = Y^exp(i0R), cl = 
•^1 — pexp(i0L) the classical Hamiltonian ^ reads 

H/N = -J^p{l -p)cose + ^{N - 1) ((1 - pf + p2) + eL(l - p) + e^p (5) 
with B = 0^-61^. 

In this paper we consider the situation of a symmetric double well where 
Cl = = Cr and initial conditions where all particles occupy the left site, i.e. |clP = 1, 
|crP = respectively p = at time t = 0. Then the corresponding conserved 
energy per particle reads E/N = H{p = 0)/N = U{N — l)/2. We assume that 
self-trapping occurs if an equal population of both sites, i.e |clP = 1/2 = |crP 
respectively p = 1/2 is no longer possible. Conservation of energy yields (cf. e.g. [1]) 
E/N = H{p^ l/2)/iV which leads to the condition cos 6* = f/(iV-l)/(2J). |cos6l| < 1 
implies that an equal occupation of both sites is only possible for \U{N — 1)/J| < 2 
so that self-trapping occurs if 

U(N-l) ^ 

A = ^—j > 2 (6) 

(for U > 0, J > 0). Please note that the effect considered here is referred to as running 
phase self-trapping and should not be confused with so-called 7r-phase self-trapping 
that is closely related to the appearance of symmetry breaking nonlinear eigenstates of 
the Gross-Pitaevskii equations ([3]), (|4]) (see e.g. [4l[T5] for a more detailed discussion). 
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3. Heuristic mean-field approach 



In the mean-field case, described by the Hamiltonian ^ and the corresponding Gross- 
Pitaevskii equations (jS]), (|4]), the site populations |clP and |crP are periodic in time 
[3]. In the following we can thus quantify the self-trapping by means of the time- 
averaged relative population imbalance (cf. [T5] ) 



2=|clP-|cr|2==1-2p -1<z<1 (7) 

where the overbar denotes time averaged quantities and p = |crP as introduced in 
the previous section. For the mean-field Gross-Pitaevskii case an exact solution in 
terms of elliptic integrals was derived in |12| . Here we instead aim at a convenient 
approximation in terms of elementary functions. 

To this end we first consider the noninteracting limit of our system, as described 
by equations ([2]), @ with U = 0. If at time t = all particles are in the left well, the 
occupation |cR(i)p = p{t) in the noninteracting two state quantum system is given by 
(see e.g. [TB] ) 



v{t) 



,P + A2 



sm 



VP + A- 
2h 



t 



(8) 



where A = eL — cr is the difference of the on-site chemical potentials. For convenience 



a brief derivation of this standard result is given in Appendix A 



In the interacting case with a finite U > the time-dependent local chemical 
potentials of the two sites read eL + U{N - l)\cL{t)\'^ and eR + U{N - l)\cR{t)\'^ 
respectively (cf. (O, ©)■ We therefore introduce the interaction into (|5]) in a heuristic 
manner by replacing A with the time-averaged difference of the on-site terms, i.e. we 
set A = EL - ER U{N - l)(|cLp - |cr|2) el - er + JA(1 - 2p). For the special 
case of a symmetric double well where el = = er equation ([5]) thus becomes 



P{t) = 



1 



1 + A2(l - 2p)2 
Self-consistency then requires 



sm 



v/l + A2(l 



1 

P=t: 



1 



2 1 + A2(l-2p)2 



(9) 



(10) 



where we have inserted the time average sin" t 
z via ([7]) , 

1 



1 



1+A22 

which yields the cubic equation 

1 



A2 



1/2. It is convenient to replace p by 



(11) 



(12) 

0. As shown in the previous section. 



For A < 2 equation (fT^ has only one solution z 
conservation of energy requires that z > (i.e. population in the left well higher than 
in the right one) if A > 2. Thus we obtain 

0, |A|<2 

|A|>2 



z{A) 



1 



(13) 
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Figure 1. Numerically exact time-averaged relative population imbalance z (*) 
as a function of the scaled interaction strength A = U{N — I)/ J compared to the 
approximation II13II (-). 



for the time-averaged relative population imbalance. The negative square root solution 
z{A) = 1/2 — (1/4 — l/A'^y/^, |A| > 2 is discarded since it docs not yield the correct 
limit z — ^ 1 for A — > oo. 

Figure [T] compares the approximation (|13p with the numerically exact values 
obtained by integrating the Gross-Pitacvskii equations ([3]), (H]). Considering the 
simplicity of our approach we observe a good agreement between the curves. The 
deviations for A > 2 are mainly due to the fact that the ansatz does not take into 
account the change in the shape of the oscillations from sinusoidal to Jacobi elliptic 
caused by the nonlinearity (cf. [3]). 

4. Heuristic semiclassical ansatz: Adding quantum fluctuations 

In this section we consider the Bosc-Hubbard dimcr ^ with a finite number N of 
particles. In [SJin] the dynamics of the initial state |iV, 0), i. e. the number state where 
all particles are in the left well, was analysed in the regime of strong interaction by 
treating the tunneling terms proportional to J as a small perturbation. It was found 
that the self-trapping effect is destroyed by quantum correlations in the long-time limit 
but that it can still be observed for a considerable time span that quickly increases 
with the number of particles in the system. Even for moderate interaction strengths, 
where the perturbative treatment of the tunneling terms is not a good approximation, 
the self-trapping effect can still be observed for a considerable time span before it is 
destroyed. This is illustrated in figure [2] which shows the time-dependent population 
imbalance z{t) for N = 100 and A = 2 for the initial state \N, 0). The time is given in 
multiples of to ~ 27r/a;o where = J/h is the Rabi frequency in the noninteracting 
case U = (cf. [Appendix A[ ). For shorter times the system shows self-trapping 
similar to the mean- field case (left panel). For longer times an oscillation with a 
higher amplitude and a longer period is revealed (right panel). For even longer times, 
which are beyond the range of numerical accuracy for the given parameters, the self- 
trapping effect is expected to be completely destroyed with the system showing a total 
population transfer between the two sites (cf. discussion in 012]). 

In the short time regime the dynamics of the system is quite well described by 
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Figure 2. Numerically exact population imbalance z(t) as a function of time for 
A'^ = 100, A = 2 and J = 1. The time is given in multiples of to = 2-Kh/J. 



semiclassical phase space methods where quantum mechanical phase representations 
of the initial quantum state like the Wigner [7] or Husimi [5J IH] distributions are 
propagated in phase space according to the Gross-Pitaevskii equations ([3]), (|4]). 

In the spirit of these semiclassical methods wc want to add the quantum 
fluctuations of the initial state to the simple mean-field description developed in 
the last section. Since the expression (jlip for the timc-avcragcd relative population 
imbalance only depends explicitly on the amplitude of the oscillation (and thus on the 
local particle numbers in the wells) but not on the relative phase between the wells we 
want to incorporate the quantum fluctuations of the initial state via the fluctuations of 
the local particle numbers. To estimate the particle number fluctuations we consider 
the approximate mean-field dynamics of the relative particle number in the right well 
given in ([5]), which can be conveniently rewritten as 

P{t) = Pamp Sin^(lUAt) (14) 

with the abbreviations pamp = (l + A2(l-2p)2)-i and wa = a/1 A2(1 - 2pYJ/{2h). 
For interaction strengths A > 2 near the transition point, where we expect the 
quantum fiuctuations to have their greatest influence, the model (fT3)) predicts z ^ 1/2 
which implies p « 1/4 and Pamp ~ 1/2 , i.e. half the particles take part in the 
oscillation. 

Particle number fluctuations of mean-field like states were recently analysed in 
[T7] . Following [T7] we define the operators 

oo = yr~paL + Vp«R > a± = -VPflL + >/! -paR (15) 

where cio destroys a state with occupations 1 — p in the left and p in the right well 
respectively and aj_ destroys a state orthogonal to the first one. A mean-field like 
state with N particles and occupations 1 — p in the left and p in the right well is thus 
given by 

\i^,) = {mr^/\alf\Q). (16) 
The inversion of relation (ITSI) reads 

The fiuctuations of the particle number fi^ — a^aR in the right well are given by 
^'"^R = (^oI^rIV'o) — (^oI^rIV'o)^- Using p7)) can be expressed in terms of clq 
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and so that the expectation values with respect to \'iPq) can be evaluated in a 
straightforward manner which yields (cf. |17j ) 

An2 =iVp(l-p). (18) 

For an initial state with p = (i.e. all particles in the left well) the particle number 
fluctuations vanish and there are only phase fluctuations which cannot be incorporated 
easily into our classical model ([T3)) whereas it is quite straightforward, as shall 
be shown in the following, to include number fluctuations. Since the approximate 
mean- field dynamics is periodic (as is the exact one, see e.g. |3]) each mean- 
field state at some time t g [0,7r/wA] creates the same dynamics (apart from an 
irrelevant phase shift) if chosen as the initial state. We therefore consider the state 
at t = tt/{2uj\) as initial state where p = Pamp and particularly Pamp = 1/2 for 
A = 2. For this state the particle number fluctuations (jl8p assume their maximum 
An|^ = Npamp{l — Pamp) = N/4:. This corresponds to a standard deviation of 
the amplitude Pamp = (^r)/^ of std(pamp) = AfiTi/N — 1/{2^/N). The standard 
deviation of the time averaged population imbalance z is then given by 

std(z) = 2std(pamp)sin2(wAt) = -4= ■ (19) 

2v N 

Applying the central limit theorem we assume the fluctuations to be Gaussian for the 
particle numbers A'^ > 50 considered in the following. 
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Now the initial quantum fluctuations can be added to our classical description by 
means of tlic replacement z — >■ z + (5z on the right hand side of 

z = l- — - , (20) 



where 5z is a Gaussian random variable with zero mean and standard deviation (|19p . 
To lowest order in 5z we obtain the cubic equation 

z^ + (2(5z - 1)5^ + ( ^ - 2(5z ) z = (21) 



which yields, in analogy to ([T5|) . 

z(A,5z)- j ,_^_^y^^T^Tj7T^^ (j,-+i)2_^>0 ■ 

The semiclassical result for the time-averaged relative population imbalance can thus 
be directly obtained by averaging over (22) with Sz ^ (/ where C is a standard 
Gaussian random variable with mean(^) = and std(C) = 1. 

Neglecting fluctuations with Sz < —1/ A — 1/2 the time-averaged relative 
population imbalance can alternatively be expressed as the integral 



m= o-^+l o+- -T^ ^=^dx (23) 




= d.,/(,+..) -^^0(f-) (24) 

where Xq = 1/A- 1/2, aN = l/(2v^), 4>ix) = cxp(-a;2/2)/\/27r is the standard 
normal distribution and $(x) = J^^dx <l){x) — [1 + erf (x/\/2)]/2 the corresponding 
cumulative distribution function. The second term in (|24p is small compared to the 
first and the third one and can be neglected. Approximating the square root term 
under the integral by its value at x* = max(xo, {xq + ctjv)/2, 0) we thus arrive at 




In figure [3] we compare our approximation for z (|25p with the numerically exact 
results for different values of the particle number N. To obtain the numerically exact 
short time results we compute the dynamics of the system in the Fock basis of states 
\N — n,n), < n < N for the initial state \N,0) and average over the time interval 
[0, lOO^o]- Generally we observe a good agreement between the approximation and 
the numerical results. The small deviations around A « 2 that are inherited from the 
heuristic mean-field approach are related to the influence of the mean-field interaction 
on the shape of the oscillations. The deviation is less pronounced for smaller particle 
numbers since the quantum fluctuations lead to deviations from the characteristic 
shape of the mean-field oscillations (cf. figure H]) that are given by Jacobi elliptic 
functions [3]. The additional deviations observed for = 50 are caused by the fact 
that for this relatively small particle number some effects of low frequency modes 
(cf. figure m and the discussion at the beginning of the section) can be felt even for 
short times for interaction strengths around A r:! 2. Qualitatively, our results confirm 
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the behaviour found in a previous numerical study |12j : The quantum fluctuations 
cause a broadening and softening of the transition region between the Josephson 
oscillation regime and the self-trapping regime. Quantitatively, however, there are 
some deviations because in |12] the time averages arc performed over much longer 
time intervals such that the system dynamics is influenced by the low frequency modes 
mentioned above. 



1.8 . ^ ^ ^ ' 
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^■^3 500 1000 

N 

Figure 4. Critical interaction strength (I26I I as a function of tlio particle number 
N for the threshold value a = 0.001. 

Since there is no longer a sharp transition point for finite particle numbers iV 
we quantify the shift of the transition point by considering the interaction strength 
A for which the population imbalance z surpasses a certain small threshold value 
a (cf. [12] )• This happens in the region where x* = xq = 1/A— 1/2, so that 
z(A) w $ ^(1 — 2/A)a//V^ /2. Setting this expression equal to a wc obtain the 
condition 

Ac, w 2 (l -iV-i/2$-i(2Q,))"^ (26) 

with the inverse function = •\/2erf~"'^(2a; — 1) known as quantile function. 

For high particle numbers S> 1 the critical interaction strength behaves like 
Aq « 2 + 2<i>^^(2Q;)A^^^/^. In the limit A^ — > oo one recovers the mean-field result 
A = 2, which is independent of a, indicating a sharp transition point. Figure |4] 
illustrates the dependence of the critical interaction strength ([^S]) on the particle 
number A^ for the threshold value a = 0.001, ^-^{2a) w -2.8782. 

In the limit A — > oo of infinite interaction strength the approximation (j25|) for 
z does not converge exactly to unity but only to a value close to one depending on 
N . This is an artifact of neglecting the fluctuations with 5z < —1/A — 1/2 which are 
only relevant for very high interaction strengths. If these fluctuations are taken into 
account equation (^5)) becomes 

and the correct limit z — >■ 1 for A — > oo is recovered. 



^ nii^)+^ hii^^) (27) 



Heuristic approach to BEC self -trapping in double wells beyond mean-field 10 
5. Summary 

In this paper we have described the transition to self-trapping of Bose-Einstein 
condensates in double wells by a single scalar quantity, namely the time-averaged 
relative population imbalance between the wells as a function of the scaled 
interaction strength. Using a heuristic ansatz we have derived convenient closed 
form approximations for the time-averaged relative population imbalance in the 
two mode Bose-Hubbard approximation that are valid in the classical (mean-field) 
and semiclassical parameter and time regimes respectively. The comparison with 
numerically exact results has revealed a good agreement. Our technically simple 
treatment complements more rigorous numerical studies of the problem based on 
semiclassical phase space methods. 
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Appendix A. The noninteracting two-mode system 

In the noninteracting case the dynamics of the site amplitudes is governed by the 
linear Schrodinger equations 

i^iCL = -^CR + eLCL (A.l) 

ihcR = --^CL + eR cr (A.2) 

which are readily obtained by setting U — in the Gross-Pitacvskii equations ([3]), 
Q. Inserting the ansatz CB,(t) = exp{iflt) and eliminating cl yields the characteristic 
equation + (el + fRj^fi + eLCR — = with the solution 



hn± = ± (A.3) 

where A = er — eL is the difference of the on-site chemical potentials. The amplitude 
in the right site can thus be written as the superposition CB,{t) = exp(ir2+f;) + 
exp(if2_i). Inserting the initial conditions ci,{t = 0) = 1, CR(t = 0) = into (|A.1[) . 
(|A.2[) we obtain ifiCY{{t = 0) = — J/2. These initial conditions for cr and cr lead to 
j4_ = — and ^4+ = J/(2a/ A^ + J^) respectively. Thus we arrive at 

cn{t) = ^J== cxp I -i^-^^t ] sin I '"l^ t I (A.4) 



V J2 + A2 

which implies the result ([5]) for p{t) — |cR(t)| 
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